Application of the Hilbert Space Average Method on Heat Conduction Models 



in 
o 
o 

(N 

> 
O 

o 

m 



J2 

o 



■4— > 

CO 

c3 



I 

C 

O 

o 



> 

(N 



in 
o 

c3 



13 

o 
o 



X 
S3 



Mathias Michel, 1 Jochen Gemmer, 2 and Giinter Mahler 1 

1 Institute of Theoretical Physics I, University of Stuttgart, Pfaffenwaldring 57, 70550 Stuttgart, German^ 
2 Physics Department, University of Osnabriick, Barbarastr. 7, 49069 Osnabriick, Germany 

(Dated: February 2, 2008) 

We analyze closed one-dimensional chains of weakly coupled many level systems, by means of 
the so-called Hilbert space average method (HAM). Subject to some concrete conditions on the 
Hamiltonian of the system, our theory predicts energy diffusion with respect to a coarse-grained 
description for almost all initial states. Close to the respective equilibrium we investigate this 
behavior in terms of heat transport and derive the heat conduction coefficient. Thus, we are able 
to show that both heat (energy) diffusive behavior as well as Fourier's law follows from and is 
compatible with a reversible Schrodinger dynamics on the complete level of description. 
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I. INTRODUCTION 

As a central problem of non-equilibrium thermody- 
namics, heat conduction has been of vital interest since 
the beginning of the 19th century. During his extensive 
investigations on heat transport in solids Fourier discov- 
ered around 1807 [l| a proportionality between the tem- 
perature gradient VT and the heat respectively energy 
current density in the material, 



j = -A VT , 



(1) 



where A is some material constant - the heat conduc- 
tivity. Concentrating on one-dimensional chains of inter- 
acting subsystems the temperature gradient can be dis- 
cretized as |VT| = AT/AL, where AT is the tempera- 
ture difference between two adjacent subsystems and AL 
their distance. Introducing the current J = Aj through 
some area A, Fourier's law can be rewritten as 



J = -n th AT 



(2) 



with Kth = AA/AL. Instead of approaching a global 
equilibrium state with some temperature T, a solid, un- 
der some external perturbation by a temperature gra- 
dient (introduced by heat baths with different tempera- 
tures) , enters a local equilibrium state only: the temper- 
ature may vary from one macroscopically small but mi- 
croscopically large part of the system to the other. The 
heat transport itself is called normal if the conductivity 
«th remains finite, whereas a diverging K t h indicates a 
ballistic transport. Having established a normal diffusive 
behavior, all further efforts on the topic of heat conduc- 
tion in the last decades have been directed towards a 
microscopic foundation of the above phenomenological 
law. 

One of the first concepts addressing the occurrence of 
regular heat conduction from a microscopic basis is the 
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Peierls-Boltzmann theory 0. Essentially, Peierls pro- 
posed a modified Boltzmann equation replacing classical 
particles by quantized quasi particles, the phonons. He 
was able to show that in order to obtain normal heat 
conduction an anharmonicity in the interaction potential 
was necessary, associated with the well known Umklapp 
processes of phonon scattering in solids. Since quantized 
normal modes are treated here as classical particles, i.e., 
as being always well localized in configuration as well as 
momentum space, this approach contains some concep- 
tual shortcomings, though. 

Another microscopic approach to heat conduction is 
the famous Green-Kubo formula, derived on the basis of 
linear response theory 0, 0, • Thereby, the current is 
taken as the response of the material to an external per- 
turbative potential, originally, the electromagnetic po- 
tential. By a direct mapping of these ideas on thermal 
transport phenomena (perturbations due to thermal gra- 
dients 0) serious questions remain open, since in this 
case it is not possible to define a proper potential in the 
Hamiltonian of the system resulting in a temperature gra- 
dient. 

To account for the problem of defining a proper poten- 
tial for a thermal perturbation, there have been efforts to 
model the system coupled to several heat baths in terms 
of an open quantum system |7j, |8|, |9j. These methods 
reveal a normal transport behavior even in very small 
quantum systems, but are numerically challenging, espe- 
cially for larger systems, since the investigations are done 
in the Liouville space instead of the Hilbert space. 

All of the above microscopic approaches to heat con- 
duction assume that the system under consideration is 
perturbed by some external environment (e.g. a scenario, 
where a solid is heated at one side cooled at the other). 
As already predicted by Einstein in his considerations on 
Brownian motion |10l 1 1 1 1] , the diffusion constant appear- 
ing in a local equilibrium steady state scenario of heat 
conduction should equal the decay constant for the re- 
laxation from a non-equilibrium state to the global equi- 
librium in the same system. In the latter case there is no 
further external perturbation involved. In the same spirit 
we will introduce yet another approach to heat transport 



2 



1 




H = 1 = 2 p = N 



FIG. 1: Heat conduction model: N coupled subunits with 
ground level and band of n equally distributed levels (AE = 
1). Black dots indicate the initial state used. 

within quantum systems. 

Recently, the relaxation behavior of small quantum 
systems has been considered within the theory of quan- 
tum thermodynamics [l2l IT3L IbH Il5j . calling for a di- 
rect solution of the Schrodinger equation. Based on the 
Hilbert space Average Method (HAM) [Tl HI [T7L fTs| . 
the complicated Schrodinger dynamics of a multi level 
system may be reduced to a set of simple rate equations. 
In the present paper we will discuss the application of 
HAM on heat conducting model systems. Furthermore, 
we will discuss the connection of heat diffusive behav- 
ior with energy diffusion and present a comparison of 
the results of HAM with the full numerical solution of 
the Schrodinger equation. The goal of the present inves- 
tigation is to show how heat respectively energy diffu- 
sive behavior may emerge from first principles - here the 
Schrodinger equation [19j . 

II. SINGLE BAND MODEL 

The class of systems we are going to analyze first is 
depicted in Fig. ^ and described by the Hamiltonian 

N N-l 

H = J2Hioc(v) + \J2 v(»>» + 1 ) ■ ( 3 ) 

fl— 1 /i— 1 

Here N identical subunits according to the same lo- 
cal Hamiltonian H\ oc (n) are assumed to have a non- 
degenerate ground state and a band of n exited states 
each, equally distributed over some band width 5e such 
that the band width is small compared to the local en- 
ergy gap AE (see Fig. ^ 5e <C AE, Se in units of AE). 
These subunits are coupled by an energy exchanging next 
neighbor interaction V(fj,, ji+ 1), chosen to be a (normal- 
ized) random hermitian matrix allowing for any possible 
transition such as to avoid any bias. Our results will turn 
out to be independent of the exact form of this matrix. 
We choose the next neighbor coupling to be weak com- 
pared to the local gap (A <C AE, A in units of AE). This 
way the full energy is approximately given by the sum of 
the local energies and those are approximately given by 
(H} oc (fi)) w AEP^, where P M is the probability to find 
the /i-th subsystem in its excited state. One of the big 
advantages of the present model is the clear partition of 
the Hamiltonian in a local and an interaction part, re- 
spectively. This partition allows for a unique definition 



of both a local energy or even a local temperature as well 
as a proper current. However, note that the presented 
method itself does not need such a special kind of model 
system but is in principle also applicable to more general 
situations. 

We now restrict ourselves to the single excitation sub- 
space, i.e., all states with one subsystem excited (irre- 
spective of which state within the band), all others in 
the ground^ state. Let us introduce a set of projection 
operators in this subspace, projecting out the set of 
states corresponding to the /i-th subsystem occupying its 
excitation band, rather than its ground state (these op- 
erators feature the unit operator of dimension n for the 
excited band of the /x-th system else zeros). Thus, the 
total system state is given by 

iv>> = E^^> = Ei^>> (4) 

where none of the is normalized individually. Fur- 
thermore, one may introduce also some standard tran- 
sition operators P^v (unit operator at position (/i,u)). 
According to the properties of transition operators we 
find P^vP^'y = ^uti'Pfiy and hence 

Pn+l,nPn,H = f^+l./x+l-F/H-l./i • (5) 

Both, the local and the interaction part of the Hamilto- 
nian (Q in our reduced Hamiltonian model can now be 
written in terms of these new operators as 

H\oc(v) = h^P^n , (6) 
V(jjl, fj, + 1) = (v^+iP^+i + vh^-Ph-i.m) ' ( 7 ) 
with a local diagonal Hamiltonian matrix and a her- 
mitian interaction matrix w^.^+i = vt+i „• 

III. TIME EVOLUTION 

The dynamics of the model system defined by the 
Hamiltonian respectively iJBJ and Q is governed by 
the Schrodinger equation. Here we will mainly be in- 
terested in a coarse-grained description suppressing the 
information about the state index within the respective 
bands, concentrating, e.g., on the total probability of oc- 
cupying a band or not, only. 

After transforming the Hamiltonian to the interaction 
picture we perform a second order Dyson expansion of 
the time evolution for the short time step r 

Mr)) « (1 - l fJ {1) ~ ^U {2) )\m) = D MO)) , (8) 
with the two time evolution operators 

U (1) (t) = [ At'V{t') , (9) 
J 

/y(2)( T ) = [ At' f T &t"V(t')V(t") . (10) 
Jo Jo 
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A pertinent quantity in the context of heat conduction 
will be the change of local energy in one subunit \x in the 
time interval r given by the probability 

P»(t) = (Mr)\Mr)) = mr)\P^P^mr)) . (11) 
Using the short time evolution operator JHJ, one gets 

P„(r) = (m\^P^P^D\m) ■■= (V>(0)|i^(0)) . 

(12) 

For simplicity we suppress in the following the explicit 
notation of the time dependence of the operators, states 
and probabilities. According to the definition of D we 
find in second order 

A — P P 4-—UWP P — — P P fl^ 
--P P U {2) - —(U i2) ) f P P 

+ ^U (1) P^P^U {1) ■ (13) 

(Note that is hermitian due to the hermiticity of the 
interaction). 

Due to the additive representation of the interaction 
Q and the definition of the first oder time evolution op- 
erator © one may find the analogous off-diagonal block 
form for the first oder time-evolution operator 

N-l 

= E ("^+1^+1 + • (14) 

The time dependent coefficients u contain the integration 
of the interaction matrix v over r' (cf. Q). 

Unfortunately, the second oder time evolution operator 
(|ll)fl has a more complicated structure in terms of the 
transition operators. Introducing the interaction Q into 
HI U| l (using the property JSJ) one finds 

= f At' f dr"J2(v^ +1 (r')v, +1 ^r")P^ 
Jo Jo 

+ f Al +i, Al (-r')«p,A'+i( T ")-PM+i,M+i 

+ U Al!Al+ l(T')u Ai +l, Al +2(T"") P ft^+2 

+ u m +i^(t')u m ,m-i( t ")4+i,m-i) ■ (I 5 ) 
This operator decomposes into a diagonal and an off- 
diagonal part = U^l g + Uf^ with 

and a suitable definition of the coefficients u^ 2 ' . By using 
(J5J one may convince oneself of the operator identities 

U (1) P^ = + P^+i)U {1) , (17) 

P^U g = U Q g! (P M+2 , M +2 + Pfj.-2.n-2) , (18) 

Piu,iJ.U di = U di P^^ . (19) 



By means of the expansion of the time-evolution opera- 
tors in terms of projectors and transition operators and 
those operator identities one is able to reformulate 
by exchanging time evolution operators and projectors as 

Ay = Pfj,,fj,Pfj,,fj, 

+ J(P () -i,,- l +P (l+1 ,, + i)P«P (t , (1 

- rPft/i^^H^-i.p-i + P/i+l.Ai+l) 

+ \P^-i{U (1) )\P^-i+P^+i) 

+ ^p M +i, /1 +i(^ (1) ) 2 (p M -i, / x-i + 

_ ip f;/( 2 ) + (f/( 2 ) ^"lp 

— tP^^U'oS {P^-241-2 + P/J.+2./1+2) 

~ ji.Plx-2.fj.-2 + P^+24i+2)(U^y Pfj.fj (20) 

Since P M = (iplA^tp) is a probability we get, furthermore, 
the normalization condition ^2 /Jj {'4 ) \A f j i \ip) = 1. This con- 
dition is already fulfilled in the zeroth order of all proba- 
bilities J2 f j('4 , \Pjl /J 1 / 1 ) — 1- D ue to this fact the sum over 
all other orders has to vanish. All first order terms to- 
gether cancel out each other directly. Remaining second 
oder terms lead to the additional conditions 

Wr-ilU™ + Wi) - (^-i|(f> (1) ) 2 |^+i) , 
(Vvl^fig + (Ugl s ) Vm> = 2(^|(L> (1) ) 2 |^) • (21) 

IV. ENERGY TRANSPORT IN THE SINGLE 
BAND MODEL 

According to l|ll|) the expectation value of the opera- 
tor (|20|l defines the time dependence of the probability 
Pp. We approximate this expectation value by an aver- 
age over some appropriate Hilbert space compartment - 
say the mean of the expectation value over all reachable 
states - called the Hilbert space average [■■•]. 

The Hilbert space average method (HAM) is, in es- 
sence, a technique to guess the value of some quantity 
defined as a function of the full system state if \ip) 
itself is not known in full detail. Here it is used to pro- 
duce a guess for the expectation values {iplA^ifb) if the 
only information about \ip) is the given set of initial ex- 
pectation values (il>\P£ \tfb) — P M (0). Such a statement 
naturally has to be a guess, since there are, in general, 
many different |V>) that are in accord with the given set of 
PfiiO) but possibly produce different values for (tplA^ip) . 
Thus, the key question for the reliability of this guess is 
whether the distribution of the A^ produced by the 
respective set of \tp)'s is broad or whether almost all those 
(tp\A fi \'ijj) are approximately equal. It can be shown that 
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if the spectral width of A^ is not too large and is high- 
dimensional almost any \tp) yields an expectation value 
close to the mean of the distribution of the (^l-A^I^) and 
thus the HAM guess will be reliable To find that 

mean one has to average with respect to the \ip)'s. This 
is called a Hilbert space average and denoted as 



P M (r) = ^(0)|4|V(0)) 



\A„ 



(22) 



Such averages are well known and derived in all detail in 
[Til lri| . For any hermitian operator A one finds 







Tr M {i} /x = 
else . 



(23) 



with the partial traces in the respective subspace /i with 
dimension (note {tp^ip^} 7^ 1). 

Plugging in the operator A^ from l|2U|l into the proba- 
bility (|12|l . approximating the expectation values by the 
Hilbert space average and furthermore using the operator 
identities (J2TJ, one finds for the probability 

^(T) = (^|^)-^I(^l(f> (1) ) 2 |^)] 

l[<^ +1 |(i>W) 2 |vv + i}] 

^Kf/W) 2 !^)]. (24) 



1 

K 2 
1 

¥ 



Furthermore, using (|23[1 with the dimension d M 
each subspace, we end up with the probability 



2 (Vvl^v 



h 2 n 

J_(Vv+i|Vv+i) 

h 2 n 
1 (^-il^-i) 
h 2 n 



Tr 



Tr 



Tr 



ti-l 



{(L>«) 2 } 
{(L7W) 2 } . 



(25) 



A first order approximation of these traces is straight- 
forward and done like in the derivation of Fermi's Golden 
Rule (cf. [15|,[X7j]) getting here 



Tr, 



{(^)) 2 } 



5c 



Defining the decay constant 



2?rA 2 n 
Me 



(26) 



(27) 



and remembering that (VvlVv) = P/j(0) = P^ one even- 
tually finds 



(28) 



Ph(t) =P tt -K(2P fi - - P^x t . 



The above approximation for the trace of the square 
of the time-evolution operator (12(:i[l is only valid within a 



small time interval like for Fermi's Golden Rule implying 
the two conditions (see |l7j l 



A 1 A 2 
T n>~, J? n«l 



(29) 



As long as these conditions are fulfilled the given second 
order time evolution approximation is valid and one may 
iterate l|28ll in the limit of r being extremely small to find 
the rate equation for the system. Since the system is a 
finite chain of N subsystems we get a slightly different 
equation for the first and the last system in the chain. 
The complete system of differential equations, thus, reads 
(2<(i<N-l) 



dPi 
dt 
dPy 
dt 
dP N 
dt 



-k(Pi-P 2 ) 
-k(p n - Pv-i 



(30) 

(31) 
(32) 



Under the conditions Ij29(l the system should follow this 
rate equation thus establishing a statistical behavior from 
pure Schrodingerian dynamics! 

From l|31|l one finds that the change of the probabil- 
ity in the /zth subsystem only depends on its left and 
right neighbor. Therefore, if the above criteria (|29|l are 
fulfilled, so that the system indeed follows the respective 
rate equation, we may concentrate on the minimal model 
of two subunits only, without loss of generality. Because 
of the special next neighbor relationship of the above dy- 
namical equation the results of the minimal model can be 
generalized to the behavior of a chain of a large number 
of such subsystems. 



V. NUMERICAL RESULTS 

That such a behavior as described in the last Sections 
is, indeed, feasible can be demonstrated by a comparison 
of the solution of the above rate equation with the exact 
solution of the full system Schrodinger equation. For this 
reason we consider a model consisting of three subunits 
with n = 500 levels each in the exited band of width 
5e = 0.05 coupled by a random energy exchanging next 
neighbor coupling A = 5 • 10~ 5 . For these parameters the 
two criteria l|29() are fulfilled. Initially we start with an 
excitation in the first system, all other systems being in 
their ground state, observing the redistribution of energy. 
The comparison for N — 3 in Fig. [21 shows a very good 
accordance between the solutions of these two different 
dynamical descriptions - both the solution of the full 
Schrodinger equation as well as the rate equation 1301) - 
(|32|l . Furthermore, Fig.[3]shows the same comparison for 
a chain of increased length N = 5. Again, one finds a 
very good accordance. 

To investigate the accuracy of the HAM for, e.g., the 
probability in the first system P\(t) we introduce D 2 , 
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FIG. 2: Probability to find the excitation in the /J, = 1, 2, 3 FIG. 4: Deviation of HAM from the exact solution (£>? < 1 

system. Comparison of the HAM prediction (lines) and the least squares) for 15000 equally distributed initial states (TV = 

exact Schrodinger solution (dots). (TV = 3, n = 500, A = 2, n = 500). 
5-HT 5 , 5e = 0.05) 



0.2 




1000 2000 3000 4000 5000 6000 
t[ft/A£] 



FIG. 3: Probability to find the excitation in the fi = 
1, 2, 3, 4, 5 system. Comparison of the HAM prediction (lines) 
and the exact Schrodinger solution (dots). (TV = 5, n = 500, 
A = 5 ■ 10~ 5 , 8e = 0.05) 




100 200 300 400 500 

n 

FIG. 5: Deviation of HAM from the exact solution (D\ least 
squares): dependence of D\ on n for TV = 3. Inset: depen- 
dence of D'f on TV for n = 500. 



as the time-averaged quadratic deviation of the exact 
(Schrodinger) result for P\(t) from the HAM prediction 

D\ = — C (Pf KM (t) ~ P°™ ct (t)fdt . (33) 
5t Jo 

To check the validity of HAM with respect to differ- 
ent pure initial states we have computed D\ for 15000 
equally distributed initial states, all corresponding to a 
two subunit system with one subunit in the upper band 
and the other subunit in its ground state. The results 
are condensed into the histogram shown in Fig.^J Obvi- 
ously, almost all initial states lead to approximately the 
same small deviation. Thus a good agreement of HAM 
with the exact result may safely be expected regardless of 
HAM being a best guess only. (In fact this is the essence 
of thermodynamical statements: They hold typically, but 
not always.) 

Furthermore, to analyze how big this "typical devi- 
ation" is, we have computed D\ for a TV = 3 system 
with the first subunit initially in an arbitrary excited 



state, for different numbers of states n in the bands. As 
shown in Fig. [5] the deviation scales like 1/n with the 
band size, i.e., vanishes in the limit of high dimensional 
subunits. This behavior does not come as a complete 
surprise since it has theoretically been conjectured and 
numerically verified in the context of equilibrium fluc- 
tuations of non-Markovian systems |13l Il5| . The inset 
of Fig. [3 shows that D\ goes down also with increasing 
number of subunits TV but then levels off. Altogether the 
new method HAM appears to be applicable even down 
to moderately sized subsystems, i.e., not too few levels n 
in the upper band. Thus, as follows from Fig. 0] a sys- 
tem with only a single excited level would lead to large 
deviations from the rate equation behavior. 

Using instead of pure initial states mixed states (which 
are typical in the context of thermodynamical phenom- 
ena) a drastic further reduction of D 2 can be expected. 
This is due to the fact that pure state fluctuations may 
cancel partially if added together. 
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VI. ENERGY DIFFUSION CONSTANT 



8e 



We concentrate in the following on a model (cf. Fig.QJ 
for which the criteria l|29[> are fulfilled. In essence, the 
used subunits should be weakly coupled to fulfill the 
given criteria, i.e., the interaction between the subunits 
has to be small compared to the local energy. The dy- 
namical behavior is then approximately governed by the 
rate equation ll3U|) - (|32|> . As argued in the last paragraph 
of Sect. llVl it is enough to consider a minimal model of 
N = 2 subunits, without loss of generality. 

This minimal system behaves according to the rate 
equation i|3U|) with the decay constant 1)27(1. The respec- 
tive equation for the second system is simply obtained by 
an exchange of the indices in l|3Ufl . 

In order to compute the energy diffusion constant in 
the given model system, we have to consider the energy 
current inside the system. The energy current is defined 
by the change of the internal energy in the two sub- 
units 



dUi 
~df 



dU 2 
dt 



(34) 



The total internal energy in such a subunit is given by 
the probability to be in the excited band of the subunit 
times the width of the energy gap, J7 M = AEP^. Thus, 
we may reformulate the current due to the change of the 
probability being in the excited band as 



■h. 



AE 



dPi 
dt 



dP 2 
dt 



(35) 



The change of the probabilities in time are given by the 
rate equation (|30[1 so that 



J h2 = -kAE(P 2 - PS 



(36) 



For a larger system with N > 2 the consideration is quite 
similar and the result in principle the same. Of course, 
the current definition is a little bit more complicated, 
since the energy relaxation in a central system leads to a 
current to both sides. But by properly splitting up the 
result into a current to the left respectively right system 
leads to the same dependence of the current between sub- 
unit [i and fi + 1 on the respective probability difference. 

According to H36|) the current is a linear function of the 
probability difference (gradient) respectively energy dif- 
ference (gradient) inside the system. This refers to a sta- 
tistical behavior since energy diffuses through the system 
only due to the energy distribution in the system. There- 
fore we may extract the energy diffusion constant from 
the above equation, identifying the parameter k with the 
transport coefficient given in (|27|) . This energy diffusion 
constant only depends on the coupling strength of the 
subunits A and on the state density of the excited band 
n/5e. (Note that in the model at hand the energy gra- 
dient is discretized as an energy difference, because of 
the discrete structure of the model, as already stated in 
Sect.Hl[cf. ©].) 



It = 2 i=2 



AE 



AE 



Ei 
n, 



:t=l®i = l 



E, 



-i = j = 



H = 2 



FIG. 6: Heat conduction model with more exited bands: each 
band consists n equally distributed levels (N = 2 subunits). 



VII. HEAT TRANSPORT IN A MULTI BAND 
MODEL 

Up to now we have restricted ourselves to the en- 
ergy transport within the system on the route from a 
non-equilibrium initial state to a global equilibrium fi- 
nal state. Of course this could be translated into heat 
transport insofar as heat is nothing but thermal energy. 
However, the states far away from equilibrium preclude a 
straight-forward interpretation: In the example of Fig. [21 
we have prepared the system initially in a state with the 
first system completely in its exited band, all others being 
in their ground states and, indeed, the system behaves 
statistically. By a simple translation to temperatures by 
fitting to a Boltzmann distribution this initial state would 
show a really strange temperature profile: one finds the 
first system at negative temperature (completely excited, 
i.e., inversion and therefore negative, infinite tempera- 
ture), all others at zero temperature. Nobody would ex- 
pect a normal heat conduction behavior in such a far 
from equilibrium situation. Nevertheless, even in this 
case we have confirmed a statistical relaxation of energy 
(see Fig.EJ. 

To investigate a situation more appropriate for heat 
transport we should use an initial state already near the 
final global equilibrium state, with only a small temper- 
ature gradient. We first start with a pure state observ- 
ing the decay of such a state by solving the Schrodinger 
equation only. Finally, one may sum up over many pure 
state trajectories finding the time evolution for an ini- 
tially mixed state. But with a temperature diffusive be- 
havior already for pure states we are sure to find it also 
for an initially mixed state, even with reduced fluctua- 
tions. 

Let us restrict ourselves in the following to only two 
subunits - the proper minimal model as argued before. 
However, let us immediately generalize the model to a 
more complicated structure with more than one excita- 
tion band per subunit. The results of HAM from the 
above considerations is valid also for this more general 
situation. The model system is shown in Fig. |SJ The 
band index i is reserved for the left system whereas j 
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belongs to the right system. We take a product state 
) = \ip\) ® | -02} as the respective initial state for all 
further investigations. For the single subsystem states 
we choose a "quasi" thermal state with temperature (3^,, 
i.e., a superposition of states according to the Boltzmann 
distribution. Here, we use the term "quasi" to label this 
special group of pure states, even if it is not a mixed 
state as usual for thermal states. Nevertheless, the occu- 
pation probability of the energy eigen states refers to the 
Boltzmann distribution! Thus, a mixture of many of such 
"quasi" thermal pure states may lead to the usual ther- 
mal state. Furthermore, we treat the bands as "quasi" 
degenerate since the band width is much smaller than 
the gap 5e <C AE. Thus, the levels inside each band are 
equally occupied according to the temperature (3^. The 
probability distribution of the first subsystem depending 
on the band energy Ei and the appropriate number of 
states n.i in the band is therefore given by 



Pi(Ei) 



77, e 



-PiEi 



(37) 




with the partition function Z(/3i). The second subsystem 
distribution is obtained simply by exchanging i by j and 
the respective quantities of the second subsystem. 

The important quantity for all further considerations 
is the probability distribution of the whole system. Since 
the total state is a product state the probability distri- 
bution of the complete system state is simply the 
product of the single probabilities 



FIG. 7: Alternative scheme to depict the complete multi- 
band spectrum: individual spectra of the two subunits on 
the Ei, Ej axes respectively; transformation to the new pa- 
rameters (E,£); black bands indicate the complete spectrum. 
Examples: (E,£) = (1,-1/2) are states, where the first sys- 
tem is in its first band, the second one in its ground state, 
(E, £) = (2, 0) are states, where both systems are in the first 
excited band etc. 



P. 



(38) 



However, it is convenient here to use a new set of coor- 
dinates for the combined spectrum in terms of the mean 
energy £ of a band and some energy displacement £ 



E := E, + E 



' 



£ ■= - (Ei 
2 V 



(39) 



(All energies, also E and £ , are measured in units of AE.) 
To support intuition we show both sets of coordinates in 
Fig. The full spectrum (black bands) can either be 
characterized in terms of the single subsystem spectra 
(E i} Ej) or by (E, £). The latter coordinates appear as a 
45°-rotation of the original ones. 

Introducing a small temperature difference AT around 
the mean temperature T so that = (&b(T ± 

AT/2)) -1 and transforming to the new parameters l|39|) 
one finds the probability distribution 



P(E,£) 



n(E,£) 
ZW)W 2 



■ exp 



ET - £AT) 
"fc B (T 2 -(AT/2) 2 ) 



(40) 



with the number of states n(E,£) in the band (E,£). 
Furthermore, in order to get a near equilibrium situation 
the temperature gradient AT has to be very small in 



comparison to the mean temperature AT <C T. There- 
fore we expand the exponential function in (|4t)(l in terms 
of AT finding in first order 



P(E,£) 



n(E,£)e k ^ 
Z 2 (T) 



1 



£ 



AT 



(41) 



where also the product of partition functions has been 
approximated in first order of the temperature gradient 
AT by Z(/3i)Z(/%) « Z(T)Z(T) = Z 2 (T). 

Again, the important quantity is the energy current 
between the two subsystems. The total energy in the 
first subunit is now given by U\ = J2iEiPi, where Pi is 
the time dependent probability distribution in the first 
system. To get the respective first system probability dis- 
tribution Pi from the time dependent full system prob- 
ability distribution we have to trace over the second 
subsystem by summing over j, finding Ui = EiPij ■ 
According to this total energy in subsystem 1 the current 
along the chain reads 



J 



d_P^ 
dt 



J dt 



or, switching to the new parameters (cf. (|39|l ) 



E,£ 



(42) 



(43) 



8 



For a single band model the sum over E would break 
down and we recover the result (f53f) . 

The change of the probability P(E, £) is again de- 
scribed by the rate equation according to the HAM. 
Since the complete model decomposes into several energy 
subspaces, the derivation is essentially just as given in 
Sect. IIVI The respective energy transport coefficient re- 
sults as a more complicated form depending on the band 
it belongs to (k = k(E, £,£')). However, we require that 
k is symmetric under the exchange of £ and £', since the 
energy transport from system 1 to system 2 should be 
as good as in the reverse direction. Therefore we end up 
with the HAM rate equation 



dP{E,£) 
di 



Y,<E,£,£')(P{E,£)-P{E,£')) (44) 



< 



16 



10 



4 6 
k B T[AE] 
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FIG. 8: Heat conductivity (1481 over temperature for a system 
with n = 500, 5e = 0.05 and A = 5 ■ 10" 5 . 



The fact that k depends on the respective band is already 
obvious from l|27|) , which depends on the state density of 
the band. Plugging in the rate equation l|44|) into the 
current l|43[> we find 

J=-^2 £k(E,£,£')(p(E,£)-P(E,£')^ . (45) 
E,e,£' 

Again, the current is simply a function of the actual state. 

For a short time step we may use the first order approx- 
imation of the "thermal" initial state 141J1 . Since we are 
summing over £ and £' the zeroth order terms according 
to H41f> will cancel. The remaining first order part 



E 

E,£,£ 



K(E,£,£')e fc B? 



k B T 2 Z 2 {T) 

(n(E,£)£ 2 - n(E, £')££') AT , (46) 

gives a heat current proportional to the temperature dif- 
ference AT. By comparing this current with Fourier's 
Law (|2J , one can extract the heat conductivity coefficient 



Kth 



E,£,£' 



n(E,£,£')e k ^ 
k B T 2 Z 2 (T) 



(n(E,£)£ 2 -n(E, £')££') 

a quantity depending only on parameters of the respec- 
tive model system and on temperature. This is a very 
general formula for the heat conductivity derived from 
HAM, valid for the class of model systems depicted in 
Fig. [S]with an arbitrary number of bands. 



VIII. HEAT TRANSPORT COEFFICIENT 

Let us return in the following to the single band model 
as shown in Fig. ^ with the partition function of a sin- 
gle n + 1-level subsystem Z(T) = 1 + n exp(-AE/k B T). 
In this case HAM leads to the unique energy diffusion 
constant k(E,£,£') = k (see $27$ ). Furthermore, in this 



single band model only the energy band E — 1AE con- 
tributes to the conductivity Kth, since for all other bands 
all possible £ parameters are zero (cf. Fig. [7J|. There 
are two possible values for the band E — AE, namely 
£ = ±AE/2 and the number of levels in these two sub- 
bands are equal n{AE, AE/2) = n(AE,-AE/2) = n, 
due to the construction of the model. One thus gets for 
the conductivity (|47H 



2irk B \ n 



nc 



AE 
k B T 



hSe 



(l + ne fc B^ ^ 



AE 



\k B T 



(48) 



The dependency of the transport coefficient Kth on 
temperature is displayed in Fig. [8] It roughly looks as 
expected for typical isolating solids. There is, however, 
no oc T 3 behavior for small temperatures, but note that 
this is only to be expected for three dimensional crystal 
lattices with a linear phonon dispersion relation for small 
momenta. The exponential behavior for small tempera- 
tures oc exp(— 1/T) results from the gap in the system 
and is difficult to treat for small temperatures since it is 
an essential singularity (typical for systems with an en- 
ergy gap). Furthermore we find a oc T~ 2 behavior rather 
than the usual oc T _1 behavior in the high temperature 
limit. Again, the usual oc T _1 is no rigorous general 
result, Peierls only states that the exponent should be 
between minus one and minus two. 

By introducing the heat capacity, the above result for 
the heat conduction coefficient can further be simplified. 
Starting with the canonical partition function of a single 
subunit Z = 1 + ne~ l3AE and the definition of the heat 
capacity in terms of the partition function 



1 d 2 In Z 
° ~ k^f 2 d(3 2 

one finds for the heat capacity of one subunit 
Mi' ' f ■ ' i \ /- N . 

k B - 



( AE \2 

I 1 + ne fc B T J 



k B T 



(49) 



(50) 
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From a comparison with i|48[l one derives 



i.e., the heat conductivity is given by the energy diffusion 
constant k multiplied by the heat capacity c of a single 
subunit. 



IX. CONCLUSION 

Rather than considering heat conduction for a system 
coupled to two heat baths with different temperature 
in a stationary local equilibrium state, we have studied 
the time evolution of a closed one-dimensional chain of 
weakly coupled many level systems. By an application of 
the new Hilbert space average method (HAM), we have 
been able to directly address a coarse-grained description. 
On this level of description the effective dynamics can be 
shown to be diffusive allowing us to derive a theoretical 
expression for the energy diffusion coefficient k from first 
principles (Schrodinger equation). Furthermore, two con- 
ditions have been established which could help to classify 
whether a given model should show a statistical diffusive 
behavior or not. 

The approximation involved (mainly the replacement 
of exact expectation values by a Hilbert space average) 
has been tested by a numerical comparison of the de- 
rived rate equation and the exact Schrodinger equation. 
Surprisingly, we have found for the present model sys- 
tem that energy diffusion occurs right from the start and 
already far from equilibrium. 

For an initial quasi thermal state featuring a small but 
finite temperature gradient we have been able to find a 
normal heat conducting behavior in terms of a finite heat 
conductivity K t h- Even for a more complicated model 
system (several bands per subunit) it was possible to de- 
rive a theoretical expression for the transport coefficient 
in terms of the parameters of the Hamiltonian. Finally, 
a reformulation of the transport coefficient has led to a 
connection with the energy transport coefficient and the 
heat capacity. 

What impact do those results for our "design model" 
have on real physical systems? For an application 
of HAM one has to organize the system as a "grain- 
structure" of weakly coupled subunits. Then, for the 



emergence of normal transport, the local bands and their 
couplings that mediate the transport should obey the cri- 
terion H'29[l individually. If this can be established HAM 
predicts that energy will diffuse from subunit to subunit. 

It has long since been appreciated that coarse-graining 
(i.e., tracing out irrelevant degrees of freedom) may 
change the effective dynamical laws even qualitatively. 
An obvious example is encountered when some acceler- 
ated electron in a solid transfers energy and momentum 
to vibrational degrees of freedom thus giving rise to par- 
ticle transport. For the origin of heat transport in insula- 
tors a correspondingly simple picture seems to be missing 
though. 

Here we have introduced a specific model class for 
which such a concept of information reduction can easily 
be incorporated. We thus have been able to clearly spec- 
ify the conditions under which the resulting effective dy- 
namics can be characterized as normal energy diffusion. 
Irreversibility and thermal resistance are thus traced back 
to a "reduced description" : "In reality" the full dynam- 
ics is still controlled by the Schrodinger equation, as has 
explicitly been checked here for still small enough total 
quantum systems. 

The desired structure may be achieved by coarse- 
graining periodic systems like spin chains, crystals, etc. 
into entities containing many elementary cells and tak- 
ing the mesoscopic entities as the subunits proper and 
the effective interactions between the entities as the cou- 
plings. Increasing the "grain size" will result in higher 
state densities within the subunits and relatively weaker 
couplings such that above a certain minimum grain size 
the criteria (|29|l will eventually be fulfilled. In this way 
a qualitatively new type of dynamical behavior can be 
seen to evolve as the observational level becomes more 
and more coarse-grained. 

Of course, the resulting subunits cannot generally be 
expected to feature the same gapped spectral structure 
as our design model. But HAM should be applicable also 
to more realistic systems, and work in that direction is 
under way. 

We thank H.-P. Breuer, M. Hartmann, M. Henrich, 
Ch. Kostoglou, H. Michel, H. Schmidt, M. Stollsteimer 
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